Zebrafish cranial neural crest

Data from: Fabian et al., Lifelong single-cell profiling of cranial neural crest diversification in zebrafish. Nat Commun 13, 13 (2022)..

Integration of scRNA-seq and snATAC-seq datasets

Load packages

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(Seurat)

## Loading required package: SeuratObject
## Loading required package: sp
## 
## Attaching package: 'SeuratObject'
## 
## The following objects are masked from 'package:base':
## 
##     intersect, t

library(Signac)
library(patchwork)
library(hexpression)

Load data

RNA

rna_datasets_paths <- list(
  dpf1.5 = "data/RNA/1.5dpf_scRNA.rds",
  dpf2 = "data/RNA/2dpf_scRNA.rds",
  dpf3 = "data/RNA/3dpf_scRNA.rds",
  dpf5 = "data/RNA/5dpf_scRNA.rds",
  dpf14 = "data/RNA/14dpf_scRNA.rds",
  dpf60 = "data/RNA/60dpf_scRNA.rds",
  dpf150 = "data/RNA/150dpf_scRNA.rds"
)

rna_datasets <- lapply(rna_datasets_paths, readRDS)

ATAC

atac_datasets_paths <- list(
  dpf1.5 = "data/ATAC/1.5dpf_snATAC_SnapATAC.rds",
  dpf2 = "data/ATAC/2dpf_snATAC_SnapATAC.rds",
  dpf3 = "data/ATAC/3dpf_snATAC_SnapATAC.rds",
  dpf5 = "data/ATAC/5dpf_snATAC_SnapATAC.rds",
  dpf14 = "data/ATAC/14dpf_snATAC_SnapATAC.rds",
  dpf60 = "data/ATAC/60dpf_snATAC_SnapATAC.rds",
  dpf210 = "data/ATAC/210dpf_snATAC_SnapATAC.rds"
)

atac_datasets <- lapply(atac_datasets_paths, readRDS)

Add annotation

The authors annotated the RNA dataset of 150 DPF (adult) zebrafish with cell types. I add the annotation to the 150dpf RNA object and then transfer it to the other datasets.

annotation <- setNames(
  c(
    "stroma 1",
    "gill progenitor 1",
    "gill cartilage",
    "gill cartilage",
    "dermal fibroblast",
    "gill cartilage",
    "gill cartilage",
    "periosteum/tendon/ligament",
    "teeth",
    "perivascular",
    "pillar",
    "smooth muscle 2",
    "gill progenitor 2",
    "smooth muscle",
    "gill stroma",
    "cycling cells",
    "tunica media",
    "stroma 2",
    "hyaline cartilage",
    "bone"
  ),
  levels(droplevels(rna_datasets$dpf150$RNA_snn_res.0.8))
)

rna_datasets$dpf150$annotation <- as.factor(unname(
  annotation[as.character(rna_datasets$dpf150$RNA_snn_res.0.8)]
))

I plot the annotation to check if it looks correct.

DimPlot(rna_datasets$dpf150, group.by = "annotation", label = T, repel = T) +
  NoAxes() +
  ggtitle("RNA annotation")

Transfer annotation to other datasets

I use the PCA transfer workflow as described in the Seurat vignette to transfer the annotation from the 150 DPF RNA dataset to the other RNA datasets. The metadata of the objects then has a column annotation with the predicted label (highest prediction score). The prediction score of that label is called prediction.score.max. Additionally, the metadata contains columns with the prediction scores for each label.

transfer_anno <- function(dataset) {
  transfer_anchors <- FindTransferAnchors(
    reference = rna_datasets$dpf150,
    query = dataset,
    features = VariableFeatures(rna_datasets$dpf150),
    reference.assay = "RNA",
    query.assay = "RNA"
  )
  predictions <- TransferData(
    anchorset = transfer_anchors,
    refdata = rna_datasets$dpf150$annotation,
    dims = 1:30
  )
  predictions <- predictions |>
    mutate(
      predicted.id = factor(
        predicted.id, levels = levels(rna_datasets$dpf150$annotation)
      )
    ) |>
    dplyr::rename(
      annotation = predicted.id
    )
  
  AddMetaData(
    dataset,
    metadata = predictions
  )
}

rna_datasets[1:6] <- lapply(rna_datasets[1:6], transfer_anno)

## Performing PCA on the provided reference using 2000 features as input.

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand

## Warning: Command RunPCA.RNA changing from SeuratCommand to SeuratCommand

## Projecting cell embeddings

## Finding neighborhoods

## Finding anchors

##  Found 197 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Performing PCA on the provided reference using 2000 features as input.

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand
## Command RunPCA.RNA changing from SeuratCommand to SeuratCommand

## Projecting cell embeddings

## Finding neighborhoods

## Finding anchors

##  Found 304 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Performing PCA on the provided reference using 2000 features as input.

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand
## Command RunPCA.RNA changing from SeuratCommand to SeuratCommand

## Projecting cell embeddings

## Finding neighborhoods

## Finding anchors

##  Found 1909 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Performing PCA on the provided reference using 2000 features as input.

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand
## Command RunPCA.RNA changing from SeuratCommand to SeuratCommand

## Projecting cell embeddings

## Finding neighborhoods

## Finding anchors

##  Found 4245 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Performing PCA on the provided reference using 2000 features as input.

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand
## Command RunPCA.RNA changing from SeuratCommand to SeuratCommand

## Projecting cell embeddings

## Finding neighborhoods

## Finding anchors

##  Found 6698 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Performing PCA on the provided reference using 2000 features as input.

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand
## Command RunPCA.RNA changing from SeuratCommand to SeuratCommand

## Projecting cell embeddings

## Finding neighborhoods

## Finding anchors

##  Found 6459 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

I plot the annotation for all RNA datasets to check if it looks correct.

(DimPlot(rna_datasets$dpf1.5, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("1.5 DPF")) +
  (DimPlot(rna_datasets$dpf2, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("2 DPF")) +
  (DimPlot(rna_datasets$dpf3, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("3 DPF")) +
  (DimPlot(rna_datasets$dpf5, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("5 DPF")) +
  (DimPlot(rna_datasets$dpf14, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("14 DPF")) +
  (DimPlot(rna_datasets$dpf60, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("60 DPF")) +
  (DimPlot(rna_datasets$dpf150, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("150 DPF")) +
  plot_layout(ncol = 2)

## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ggsave("figures/RNA_annotation.pdf")

## Saving 7.25 x 12 in image

## Warning: ggrepel: 3 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

The early stages are dominated by cycling cells. We should check if this is expected. There is not much information about this in this paper, the earlier samples are just annotated as mesenchyme. In later stages, more defined cell types appear and the annotation looks reasonable.

We should still take some time to check with literature marker genes.

Transfer annotation and clusters to RNA to ATAC

Now, I will transfer the annotation and clusters from the RNA datasets to the ATAC datasets. I will use the CCA transfer workflow as described in the Seurat vignette. I first try it with the clusters of the 1.5 DPF RNA dataset, to test if it works.

transfer_anchors <- FindTransferAnchors(
  reference = rna_datasets$dpf1.5,
  query = atac_datasets$dpf1.5,
  features = VariableFeatures(rna_datasets$dpf1.5),
  reference.assay = "RNA",
  query.assay = "ACTIVITY",
  reduction = "cca"
)

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 4529 anchors

cluster_predictions <- TransferData(
  anchorset = transfer_anchors,
  refdata = rna_datasets$dpf1.5$RNA_snn_res.0.8,
  weight.reduction = atac_datasets$dpf1.5[["lsi"]],
  dims = 2:30
)

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

dpf1.5_test <- AddMetaData(
  atac_datasets$dpf1.5,
  metadata = cluster_predictions
)
dpf1.5_test$predicted.id <- factor(
  dpf1.5_test$predicted.id,
  levels = c("0", "1", "2", "4", "6", "7", "9", "10")
)

I plot the clusters and the expression/gene activity score of some marker genes that are specific for these clusters to check if the transfer result looks reasonable.

(
  DimPlot(rna_datasets$dpf1.5, group.by = "RNA_snn_res.0.8", label = T) +
    ggtitle("RNA") |
    DimPlot(dpf1.5_test, group.by = "predicted.id", label = T) +
    ggtitle("ATAC")
) /
  (
    HexPlot(rna_datasets$dpf1.5, "pitx2", bins = 100) |
      HexPlot(dpf1.5_test, "pitx2", bins = 100)
  ) /
  (
    HexPlot(rna_datasets$dpf1.5, "hoxd4a", bins = 100) |
      HexPlot(dpf1.5_test, "hoxd4a", bins = 100)
  ) /
  (
    HexPlot(rna_datasets$dpf1.5, "hand2", bins = 100) |
      HexPlot(dpf1.5_test, "hand2", bins = 100)
  ) /
  (
    HexPlot(rna_datasets$dpf1.5, "pou3f3a", bins = 100) |
      HexPlot(dpf1.5_test, "pou3f3a", bins = 100)
  )

ggsave("figures/RNA_ATAC_integration_1.5dpf_expression_activity.pdf")

## Saving 6.5 x 12 in image

This looks reasonable. I now transfer the annotation and the clusters pairwise from each RNA dataset to the corresponding ATAC dataset.

Similar to the annotation transfer in the RNA, the metadata of the objects then has a column annotation with the predicted label (highest prediction score) and a column RNA_cluster with the predicted cluster from the RNA data. The prediction scores of these labels are called annotation_prediction.score.max and RNA_cluster_prediction.score.max, respectively. Additionally, the metadata contains columns with the prediction scores for each label.

transfer_anno_to_ATAC <- function(dataset_atac, dataset_rna) {
  transfer_anchors <- FindTransferAnchors(
    reference = dataset_rna,
    query = dataset_atac,
    features = VariableFeatures(dataset_rna),
    reference.assay = "RNA",
    query.assay = "ACTIVITY",
    reduction = "cca"
  )
  cluster_predictions <- TransferData(
    anchorset = transfer_anchors,
    refdata = dataset_rna$RNA_snn_res.0.8,
    weight.reduction = dataset_atac[["lsi"]],
    dims = 2:30
  ) |>
    mutate(
      predicted.id = factor(
        predicted.id, levels = levels(dataset_rna$RNA_snn_res.0.8)
      )
    ) |>
    dplyr::rename(
      RNA_cluster = predicted.id
    ) |>
    dplyr::rename_with(
      .fn = function(x) sub("prediction", "RNA_cluster_prediction", x)
    )
  annotation_predictions <- TransferData(
    anchorset = transfer_anchors,
    refdata = dataset_rna$annotation,
    weight.reduction = dataset_atac[["lsi"]],
    dims = 2:30
  ) |>
    mutate(
      predicted.id = factor(
        predicted.id, levels = levels(dataset_rna$annotation)
      )
    ) |>
    dplyr::rename(
      annotation = predicted.id
    ) |>
    dplyr::rename_with(
      .fn = function(x) sub("prediction", "annotation_prediction", x)
    )
  dataset_atac <- AddMetaData(
    dataset_atac,
    metadata = cluster_predictions
  )
  dataset_atac <- AddMetaData(
    dataset_atac,
    metadata = annotation_predictions
  )
  dataset_atac
}

atac_datasets <- Map(transfer_anno_to_ATAC, atac_datasets, rna_datasets)

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 4529 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 6963 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 7208 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 6387 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 9741 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 3782 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Warning: Command ScaleData.RNA changing from SeuratCommand to SeuratCommand

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 4410 anchors

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

I plot the results for all datasets to check if the transfer worked well.

(DimPlot(rna_datasets$dpf1.5, group.by = "RNA_snn_res.0.8", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 1.5 DPF")) +
  (DimPlot(atac_datasets$dpf1.5, group.by = "RNA_cluster", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 1.5 DPF")) +
  (DimPlot(rna_datasets$dpf2, group.by = "RNA_snn_res.0.8", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 2 DPF")) +
  (DimPlot(atac_datasets$dpf2, group.by = "RNA_cluster", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 2 DPF")) +
  (DimPlot(rna_datasets$dpf3, group.by = "RNA_snn_res.0.8", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 3 DPF")) +
  (DimPlot(atac_datasets$dpf3, group.by = "RNA_cluster", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 3 DPF")) +
  (DimPlot(rna_datasets$dpf5, group.by = "RNA_snn_res.0.8", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 5 DPF")) +
  (DimPlot(atac_datasets$dpf5, group.by = "RNA_cluster", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 5 DPF")) +
  (DimPlot(rna_datasets$dpf14, group.by = "RNA_snn_res.0.8", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 14 DPF")) +
  (DimPlot(atac_datasets$dpf14, group.by = "RNA_cluster", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 14 DPF")) +
  (DimPlot(rna_datasets$dpf60, group.by = "RNA_snn_res.0.8", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 60 DPF")) +
  (DimPlot(atac_datasets$dpf60, group.by = "RNA_cluster", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 60 DPF")) +
  (DimPlot(rna_datasets$dpf150, group.by = "RNA_snn_res.0.8", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 150 DPF")) +
  (DimPlot(atac_datasets$dpf210, group.by = "RNA_cluster", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 210 DPF")) +
  plot_layout(ncol = 2)

ggsave("figures/RNA_ATAC_integration_clusters.pdf")

## Saving 7.25 x 18 in image

(DimPlot(rna_datasets$dpf1.5, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 1.5 DPF")) +
  (DimPlot(atac_datasets$dpf1.5, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 1.5 DPF")) +
  (DimPlot(rna_datasets$dpf2, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 2 DPF")) +
  (DimPlot(atac_datasets$dpf2, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 2 DPF")) +
  (DimPlot(rna_datasets$dpf3, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 3 DPF")) +
  (DimPlot(atac_datasets$dpf3, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 3 DPF")) +
  (DimPlot(rna_datasets$dpf5, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 5 DPF")) +
  (DimPlot(atac_datasets$dpf5, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 5 DPF")) +
  (DimPlot(rna_datasets$dpf14, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 14 DPF")) +
  (DimPlot(atac_datasets$dpf14, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 14 DPF")) +
  (DimPlot(rna_datasets$dpf60, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 60 DPF")) +
  (DimPlot(atac_datasets$dpf60, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 60 DPF")) +
  (DimPlot(rna_datasets$dpf150, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("RNA 150 DPF")) +
  (DimPlot(atac_datasets$dpf210, group.by = "annotation", label = T, repel = T) + NoAxes() + NoLegend() + ggtitle("ATAC 210 DPF")) +
  plot_layout(ncol = 2)

## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 3 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ggsave("figures/RNA_ATAC_integration_annotation.pdf")

## Saving 7.25 x 18 in image

## Warning: ggrepel: 4 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

Save the results

Seurat objects

I save the Seurat objects as lists of objects, for the RNA and ATAC datasets, respectively.

saveRDS(rna_datasets, "saved_data/rna_datasets_list.rds")
saveRDS(atac_datasets, "saved_data/atac_datasets_list.rds")

Metadata

I export the metadata of the RNA and ATAC datasets as TSV files, with the cell names as first column. The metadata contains the (transferred) annotation and the original clusters (RNA_snn_res.0.8) for the RNA datasets, and the transferred annotation and RNA clusters for the ATAC datasets.

walk(names(rna_datasets), function(name) {
  write.table(
    rna_datasets[[name]]@meta.data |>
      rownames_to_column(var = "cell"),
    file = paste0("saved_data/rna_", name, "_metadata.tsv"),
    sep = "\t",
    quote = F,
    row.names = F
  )
})

walk(names(atac_datasets), function(name) {
  write.table(
    atac_datasets[[name]]@meta.data |>
      rownames_to_column(var = "cell"),
    file = paste0("saved_data/atac_", name, "_metadata.tsv"),
    sep = "\t",
    quote = F,
    row.names = F
  )
})

Session info

sessionInfo()

## R version 4.3.3 (2024-02-29)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 24.04.2 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /work/ntrost/miniconda3/envs/r4Seurat5/lib/libopenblasp-r0.3.28.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] hexpression_1.0.1  patchwork_1.3.0    Signac_1.14.0      Seurat_5.1.0      
##  [5] SeuratObject_5.0.2 sp_2.1-4           lubridate_1.9.4    forcats_1.0.0     
##  [9] stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5       
## [13] tidyr_1.3.1        tibble_3.2.1       ggplot2_3.5.1      tidyverse_2.0.0   
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3      rstudioapi_0.17.1       jsonlite_1.8.9         
##   [4] magrittr_2.0.3          spatstat.utils_3.1-1    farver_2.1.2           
##   [7] rmarkdown_2.29          ragg_1.3.3              zlibbioc_1.48.0        
##  [10] vctrs_0.6.5             ROCR_1.0-11             Rsamtools_2.18.0       
##  [13] spatstat.explore_3.3-3  RCurl_1.98-1.16         RcppRoll_0.3.1         
##  [16] htmltools_0.5.8.1       sass_0.4.9              sctransform_0.4.1      
##  [19] parallelly_1.40.1       KernSmooth_2.23-24      bslib_0.8.0            
##  [22] htmlwidgets_1.6.4       ica_1.0-3               plyr_1.8.9             
##  [25] plotly_4.10.4           zoo_1.8-12              cachem_1.1.0           
##  [28] igraph_2.0.3            mime_0.12               lifecycle_1.0.4        
##  [31] pkgconfig_2.0.3         Matrix_1.6-5            R6_2.5.1               
##  [34] fastmap_1.2.0           GenomeInfoDbData_1.2.11 fitdistrplus_1.2-1     
##  [37] future_1.34.0           shiny_1.10.0            digest_0.6.37          
##  [40] colorspace_2.1-1        S4Vectors_0.40.2        tensor_1.5             
##  [43] RSpectra_0.16-2         irlba_2.3.5.1           textshaping_0.3.7      
##  [46] GenomicRanges_1.54.1    labeling_0.4.3          progressr_0.15.1       
##  [49] fansi_1.0.6             spatstat.sparse_3.1-0   timechange_0.3.0       
##  [52] httr_1.4.7              polyclip_1.10-7         abind_1.4-5            
##  [55] compiler_4.3.3          withr_3.0.2             BiocParallel_1.36.0    
##  [58] viridis_0.6.5           fastDummies_1.7.4       hexbin_1.28.5          
##  [61] MASS_7.3-60.0.1         tools_4.3.3             lmtest_0.9-40          
##  [64] httpuv_1.6.15           future.apply_1.11.2     goftest_1.2-3          
##  [67] glue_1.8.0              nlme_3.1-165            promises_1.3.2         
##  [70] grid_4.3.3              Rtsne_0.17              cluster_2.1.8          
##  [73] reshape2_1.4.4          generics_0.1.3          gtable_0.3.6           
##  [76] spatstat.data_3.1-4     tzdb_0.4.0              data.table_1.15.2      
##  [79] hms_1.1.3               XVector_0.42.0          utf8_1.2.4             
##  [82] BiocGenerics_0.48.1     spatstat.geom_3.3-4     RcppAnnoy_0.0.22       
##  [85] ggrepel_0.9.6           RANN_2.6.2              pillar_1.9.0           
##  [88] spam_2.11-0             RcppHNSW_0.6.0          later_1.4.1            
##  [91] splines_4.3.3           lattice_0.22-6          survival_3.7-0         
##  [94] deldir_2.0-4            tidyselect_1.2.1        Biostrings_2.70.1      
##  [97] miniUI_0.1.1.1          pbapply_1.7-2           knitr_1.49             
## [100] gridExtra_2.3           IRanges_2.36.0          scattermore_1.2        
## [103] stats4_4.3.3            xfun_0.49               matrixStats_1.4.1      
## [106] stringi_1.8.4           lazyeval_0.2.2          yaml_2.3.10            
## [109] evaluate_1.0.1          codetools_0.2-20        cli_3.6.3              
## [112] uwot_0.1.16             systemfonts_1.1.0       xtable_1.8-4           
## [115] reticulate_1.40.0       munsell_0.5.1           jquerylib_0.1.4        
## [118] Rcpp_1.0.13-1           GenomeInfoDb_1.38.1     globals_0.16.3         
## [121] spatstat.random_3.3-2   png_0.1-8               spatstat.univar_3.1-1  
## [124] parallel_4.3.3          dotCall64_1.2           bitops_1.0-9           
## [127] listenv_0.9.1           viridisLite_0.4.2       scales_1.3.0           
## [130] ggridges_0.5.6          crayon_1.5.3            leiden_0.4.3.1         
## [133] rlang_1.1.4             fastmatch_1.1-6         cowplot_1.1.3